if !d.name eq 'PS' then begin
    device,xsize=18,ysize=8,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end
!x.margin=[8,2]
!p.charsize=1.5

if n_elements(bx) eq 0 then begin

 tsu2 = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
 tsv2 = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
 tsw2 = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)
 tsp = (total((total(p,1)/m),1)/l)
 tsth = (total((total(ptemp,1)/m),1)/l)

 nn = n_elements(tsu2)
 urms = sqrt(tsu2+tsv2+tsw2)
 year = 365.*24.*60.*60.
 time = (nxaver*dt*indgen(nn))/year

 plot,time,urms,xtitle='!8t!6 (years)',ytitle='!8u!D!6rms!N',thick=3
 plot,time,tsp,xtitle='!8t!6 (years)',ytitle='!8p!6',thick=3
 plot,time,tsth,xtitle='!8t!6 (years)',ytitle='!4h!6',thick=3

endif else begin

 tsu2 = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
 tsv2 = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
 tsw2 = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)
 tsbx2 = (total((total(bx2,1)/m),1)/l)-(total((total(bx^2,1)/m),1)/l)
 tsby2 = (total((total(by2,1)/m),1)/l)-(total((total(by^2,1)/m),1)/l)
 tsbz2 = (total((total(bz2,1)/m),1)/l)-(total((total(bz^2,1)/m),1)/l)
 tsp = (total((total(p,1)/m),1)/l)
 tsth = (total((total(ptemp,1)/m),1)/l)

 uu2=(total((total((w^2+v^2+u^2),1)/m),1)/l)
 bb2=(total((total((bx^2+by^2+bz^2),1)/m),1)/l)

 nn = n_elements(tsu2)
 urms = sqrt(tsu2+tsv2+tsw2)
 brms = sqrt(tsbx2+tsby2+tsbz2)
 year = 365.*24.*60.*60.
 time = (nxaver*dt*indgen(nn))/year
 bbx=total(bx[18,*,*],2)/l
; computing energies
 Eeq=0.5*mean(urms[300:*])^2
 Eu=0.5*uu2/Eeq
 Eb=bb2/(70.*4e-7*!pi)/Eeq

 plot,time,bbx^2/(70.*4e-7*!pi)/Eeq,xtitle='!8t!6 [yr]',ytitle='!8<E!DB!D!7u!8!N!N>/E!Deq!N!6  !N',thick=6, $
              xstyle=1,ystyle=1,yr=[0,0.1],xrange=[20,52]

endelse 


!p.multi=0
end
